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Abstract.  Results  are  presented  from  a parallel  DNS  combustion  code 
called  SENGA.  The  code  solves  a fully  compressible  reacting  flow  in  three 
dimensions.  High  accuracy  numerical  schemes  have  been  employed  which 
are  explicit  10th  order  central  finite  differences  in  space,  a third  order  ex- 
plicit Runge-Kutta  method  in  time  and  parallel  implementation  is  achieved 
through  the  Message  Passing  Interface  (MPI).  Turbulence  is  generated  nu- 
merically for  a 1283  simulation  with  Re  = 30  and  a 3843  simulation  with 
Re  = 130.  Finally,  results  are  presented  and  discussed  for  simulations  with 
different  initial  non-dimensional  turbulence  intensities  ranging  from  5 to 
23. 


1.  Introduction 

A difficulty  that  arises  in  practical  turbulent  combustion  processes,  such  as 
the  combustor  sections  of  jet  engines  and  internal  combustion  engines,  is 
the  strong  coupling  between  turbulence,  chemical  kinetics  and  heat  release. 
These  interactions  are  generally  three  dimensional  and  time  dependent,  and 
are  not  easily  accessible  to  experimental  investigation.  Therefore,  to  extract 
valuable  information  from  these  and  other  turbulent  flows,  highly  accurate 
numerical  solutions  are  required,  such  as  those  obtained  using  Direct  Nu- 
merical Simulation  (DNS). 

DNS  is  now  a useful  and  well  established  research  tool  in  the  field  of  tur- 
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bulent  combustion.  The  solution  of  the  full  governing  equations  without 
the  need  for  any  form  of  turbulence  modelling  enables  the  collection  of 
a great  deal  of  information  on  turbulent  flame  structure  and  propagation 
that  is  obtainable  by  no  other  means.  However,  to  obtain  representative 
simulations,  the  DNS  must  be  able  to  resolve  the  smallest  length  and  time 
scales  applicable  to  the  problem  of  interest.  High  resolution  and  the  need  to 
solve  in  three  spatial  directions  and  time  comes  at  a cost,  therefore  parallel 
supercomputers  are  a pre-requisite  to  carry  out  these  simulations.  Due  to 
the  high  computational  costs,  it  is  important  to  use  high  accuracy  spatial 
discretisation  schemes.  The  most  accurate  scheme  is  the  Fourier  spectral 
method,  however,  these  schemes  are  restricted  to  problems  with  periodic 
boundary  conditions  which  in  combustion  DNS  is  a restriction  since  reac- 
tants must  be  able  to  enter  and  products  and  heat  must  be  able  to  leave 
the  computational  domain.  Therefore  alternative  schemes  with  spectral  like 
accuracy  and  non  restrictive  boundary  conditions  are  required,  the  most 
popular  scheme  for  combustion  being  the  compact  scheme  (Lele,  1992).  In 
the  present  work,  10th  order  explicit  finite  differences  are  employed.  These 
offer  near  spectral  accuracy  along  with  inflowing  and  outflowing  boundary 
conditions  and  are  more  efficient  for  parallel  implementation.  (Prosser  and 
Cant,  1998)  compared  various  spatial  schemes  applicable  to  combustion 
DNS,  and  showed  the  explicit  scheme  to  perform  well  in  terms  of  accuracy 
when  the  resolution  is  known  to  be  sufficient,  and  in  terms  of  CPU  time 
and  storage  the  explicit  schemes  proved  least  expensive  on  both  counts. 
Practical  combustion  problems  have  chemical  and  turbulent  scales  and  in 
most  cases  it  is  the  smallest  chemical  scales  within  the  reaction  zone  of  the 
flame  that  are  of  interest,  not  the  smallest  turbulent  scales.  This  is  true 
for  example  in  the  spark  ignition  engine  where  combustion  takes  place  in 
a reasonably  homogeneous  mixture  of  fuel  vapour  and  air,  and  where  the 
turbulence  is  of  moderate  intensity.  These  kind  of  systems  operate  within 
the  laminar  flamelet  regime  of  turbulent  premixed  combustion  (Libby  and 
Williams,  1994).  Here  the  thickness  of  the  flame  sheet  remains  smaller  than 
the  Kolomogorov  scale  of  turbulence  and  the  Kilmov-Williams  criteria  is 
satisfied  (Williams,  1985).  The  flame  sheet  retains  the  structure  of  a lam- 
inar flame,  even  though  it  is  wrinkled  by  the  surrounding  turbulence.  The 
flame-turbulence  interactions  in  the  present  work  are  in  the  above  regime 
and  therefore  the  flame  structure  and  not  the  turbulence  defines  the  resolu- 
tion requirement.  Research  in  combustion  DNS  as  come  a long  way  in  these 
regimes  (Poinsot  et  al.,  1996),  especially  for  planar  flames.  The  present 
work  uses  flame  kernels  in  a turbulent  environment  which  have  a practical 
relation  to  ignition  problems.  The  growth  of  the  flame  is  laminar  once  es- 
tablished in  the  early  stages,  and  then  propagates  spherically  outwards.  Its 
motion  is  accelerated  by  flow  divergence  due  to  thermal  expansion  and  it 
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begins  to  interact  with  the  surrounding  turbulence  becoming  wrinkled  and 
increasing  in  surface  area.  Flame  kernels  in  conjuction  with  non-reflecting 
outflow  boundaries  on  all  faces  using  the  NSCBC  formalism  (Poinsot  and 
Lele,  1992),  allow  the  flame  structure  to  be  observed  away  from  the  influ- 
ence of  any  boundary  impositions. 

The  aim  of  the  present  research  is  to  use  a parallel  DNS  code  called  SENGA, 
developed  in  Cambridge,  to  study  flame  kernels  in  a turbulent  environment. 
The  DNS  code  solves  the  fully  compressible  reacting  flow  equations  in  three 
dimensions  and  time  with  heat  release.  The  paper  is  organised  as  follows. 
The  equations  governing  the  flow  are  outlined  in  non-dimensional  format, 
the  numerical  procedure  is  explained  and  finally  the  simulation  method  and 
results  are  presented  for  various  flame  turbulent  interactions. 

2.  Governing  Equations 

The  governing  equations  that  describe  the  motion  of  a reacting  gas  are  the 
three  dimensional  equations  for  mass,  momentum,  energy  and  a species  con- 
servation equation  (Williams,  1985).  In  non-dimensional  form  these  equa- 
tions become. 
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where  x and  t are  the  space  and  time  coordinates,  p is  the  density,  u is 
the  velocity  , P is  the  pressure,  E is  a stagnation  internal  energy,  c is  a 
reaction  progress  variable,  formulated  as  a normalised  mass  fraction  and 
rising  monotonically  from  zero  in  the  unburned  state  to  unity  in  the  fully 
burned  products  and  rki  is  the  viscous  stress  tensor 

The  non-dimensional  parameters  in  the  above  equations  are  the  Mach  num- 
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ber  M,  the  Reynolds  Number  Re,  the  Prandtl  number  Pr  and  the  Schmidt 
number  Sc.  For  a full  description  of  the  non-dimensionisation  and  the  extra 
terms  needed  to  close  the  problem,  see  (Jenkins  and  Cant,  1999). 


3.  Numerical  Procedure 

Turbulent  combustion  DNS  requires  spatial  schemes  with  high  accuracy. 
Therefore,  all  first  and  second  spatial  derivatives  are  discretised  using  a 
central  10tfl  order  explicit  scheme  obtained  from  the  central  finite  differ- 
ence approximation 


m/2 

fi  = 2 ~ fo-i) 


where  m is  the  order  of  the  approximation  (always  even).  Values  of  con- 
stants dj  are  obtained  by  expanding  in  Taylor  series  and  equating  coeffi- 
cients of  successive  orders  in  h. 

These  finite  difference  approximations  have  a stencil  width  of  eleven  points 
which  implies  5 points  are  required  at  the  boundary.  These  boundary  points 
are  again  treated  with  explicit  finite  differences  of  decreasing  accuracy  to 
the  boundary.  This  is  covered  in  detail  by  (Jenkins  and  Cant,  1999).  Time 
stepping  is  carried  out  by  an  explicit  third  order  Runge-Kutta  method 
(Wray,  1990).  This  method  requires  three  sub  steps  for  each  main  time  ad- 
vancement step  and  requires  only  two  storage  locations,  one  for  the  time 
derivative  and  one  for  the  dependent  variable. 

An  initial  field  of  isotropic  turbulence  is  generated  numerically  to  satisfy 
the  continuity  constraint  for  incompressible  flow.  The  general  requirement 
for  this  turbulent  field  is  the  specification  of  an  energy  spectrum  E(k), 
where  k is  a wavenumber  in  the  Fourier-space  representation  of  the  turbu- 
lent velocity  field.  Various  energy  specta  exist  for  this  work  such  as  (Lee  and 
Reynolds,  1985)  which  as  previously  been  used  by  the  authors.  However, 
the  present  work  uses  the  spectrum  of  (Schumann  and  Patterson,  1978)  af- 
ter studies  showed  a slower  decay  in  turbulent  kinetic  energy  compared  to 
Lee  and  Reynolds.  This  is  of  interest,  especially  for  higher  turbulence  inten- 
sity cases,  were  the  flame  interacts  longer  at  the  higher  intensities.  Finally, 
boundary  conditions  are  implemented  using  the  Navier  Stokes  Character- 
istic Boundary  Conditions  (NSCBC)  (Poinsot  and  Lele,  1992)  and  parallel 
implementation  is  achieved  through  the  Message  Passing  Interface  (MPI). 
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4.  Simulation  Results 

Flame  kernel  interactions  for  various  turbulent  intensities  will  be  presented. 
Two  simulations  were  undertaken,  one  with  a grid  size  3843  and  the  other 
with  a grid  size  1283.  Both  simulations  were  undertaken  using  64  proces- 
sors on  the  Hitachi  SR2201  at  Cambridge  University  and  table  2 shows  the 
problem  parameters  for  each  case. 


TABLE  1.  Simulation  Parameters 


Case 

Re 

Pr 

Sc 

M 

0 

a 

r 

B 

1283 

30 

0.7 

0.7 

0.00142 

6.0 

5.0 

0.83 

1225.0 

3843 

130 

0.7 

0.7 

0.00142 

6.0 

5.0 

0.83 

1225.0 

In  both  cases,  the  simulations  were  initialised  by  a Gaussian  distribution 
of  the  progress  variable  c decomposed  on  the  64  processors  with  an  initial 
turbulent  field.  All  boundary  conditions  were  non-reflecting  outflow  types 
with  a time  step  At  = 1.0  x 10~6  which  is  restricted  by  the  acoustic  Courant 
stability  criteria.  Initial  studies  on  small  grids  (643)  have  been  undertaken 
with  M = 0.0142,  enabling  larger  time  steps  At  = 1.0  x 10~4.  Compu- 
tational performance  is  reasonable,  with  the  initial  turbulent  fields  taking 
approximately  8 and  20  minutes  to  generate  for  1283  and  3843  simulations 
respectively,  and  time  stepping  takes  approximately  15  and  50  seconds  for 
each  case.  A full  parallel  performance  test  can  be  found  in  (Jenkins  and 
Cant,  1999). 

Figures  la  - lc  show  results  from  initial  values  of  5,  10  and  23  re- 
spectively. In  each  figure,  three  sets  of  results  are  plotted  at  2,  4,  6 and  8 
thousand  time  steps.  The  left  and  centre  columns  show  a 2D  slice  from  a 
laminar  and  turbulent  flame  and  the  right  column  shows  the  isosurface  of 
progress  variable  c,  all  at  the  same  time  levels  respectively.  The  laminar  and 
turbulent  flame  cases  are  plotted  with  progress  variable  contours  ranging 
from  1 (fully  burned)  at  the  centre  to  0 (unburned)  on  the  outer  contour. 
In  each  case  the  flames  can  clearly  be  seen  to  propagate  outwards.  As  the 
flames  develop  the  effect  of  turbulence  on  them  becomes  less  significant, 
since  the  turbulence  is  decaying.  From  these  figures,  an  important  feature 
of  the  flame  is  clear.  Small  scale  wrinkles  in  the  flame,  more  noticeable 
at  ^ =23,  do  not  remain.  They  appear  to  be  smoothed  out  as  the  flame 
propagates.  Also  evident  at  high  turbulence  intensity  is  the  flame  breaking 
up  and  re-joining  as  the  turbulence  decays. 

Figures  2a  and  2b  show  a (1283)  flame  evolution  at  the  centre  x-y  plane 
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(k=64)  and  the  quarter  x-y  plane  (k=96)  over  time  starting  at  the  top 
left  and  going  to  the  bottom  right  of  each  figure.  Two  interesting  features 
occur  in  these  figures,  the  first  being  an  increase  in  burned  area  with  time 
about  the  initial  ignition  point  and  the  second  is  the  evidence  of  holes  in 
the  flame  due  to  strain.  The  holes  are  more  noticeable  in  figure  2b.  Here  the 
image  plane  is  initially  at  the  edge  of  the  flame  and  as  the  flame  progresses 
holes  and  break  up  are  evident.  Figure  3 looks  at  an  enlarged  portion  of 
the  centre  plane.  Here  the  flame  front  is  moving  to  the  right  and  upwards. 
Velocity  vectors  are  superimposed  on  to  the  flame  contours  to  show  the  ef- 
fect of  flame  curvature  due  to  the  eddies  in  the  turbulent  field.  Also  at  the 
points  of  maximum  curvature  the  flame  thickness  seems  to  be  increasing 
due  to  a higher  burning  rate,  since  the  surface  area  is  increasing. 

Finally,  figure  4 shows  an  isosurface  of  the  progress  variable  at  c = 0.8.  The 
image  shown  is  for  the  inner  8 processors  of  a 64  processor  simulation.  This 
is  evident  by  the  rings  on  each  surface  due  to  the  flame  leaving  one  processor 
and  joining  the  next.  This  is  a good  indication  of  the  veladity  of  the  parallel 
implementation  as  well  as  visualising  in  three  dimensions  the  flame  surface. 


5.  Figures 


(a)  £=5 


(b)  £=10 


=23 


Figure  1.  laminar,  turbulent  and  isosurfaces  for  3 turbulent  cases  t = 2k  - 8k 
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(a)  Centre  Plane  (b)  Quarter  Plane 


Figure  2.  2D  slice  at  centre  x-y  plane  and  quarter  x-y  plane 


Figure  3.  Enlarged  2D  slice  of  progress  variable  contours  and  velocity  vectors 


Figure  4 


Isosurface  of  progress  variable,  c = 0.8 
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